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Abstract 

A Kalman filter can be used to determine material parameters using uncertain ex¬ 
perimental data. However, starting with inappropriate initial values for material 
parameters might include false local attractors or even divergence. Also, inappro¬ 
priate choices of covariance errors of initial state, present state, and measurements 
might affect the stability of the prediction. The present method suggests a simple 
way to predict the parameters and the errors, required to start Kalman filter based 
on known parameters that are used to generate the data with different noises used 
as “measurement data”. The method consists of two steps. First, an appropriate 
range of parameter values is chosen based on a graphical representation of the 
mean square error. Second, the Kalman filter is used based on the selected range 
and the suggested parameters and errors. The method of the filter significantly 
reduces the iteration time, and covers a wide range of initial suggested values for 
the parameters compared with the standard Kalman filter. When the methodology 
is applied to real data, very good results are obtained. Diffusion coefficient for 
bovine bone is chosen to be a case study in this work. 
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1. Introduction 


The Kalman filter is an inverse method to determine variables or parameters 
using input data with more noise and get output data with less noise. It is firstly 
presented by R.E. Kalman [[11 in 1960. Kalman filter has the advantages of tak¬ 
ing the random noise for state and measurements into eonsideration, also it is an 
optimal estimator for linear models beeause it minimizes the mean square error 
between the state. In addition, it eonverges quiekly. A more eomplete introduetion 
to the Kalman filter is given by Brown [[3. The Kalman filter ean be found under 
different updated forms that used in many different fields sueh as traeking objeets 
[[SHSI, eontrol systems [l6l|71, and weather foreeast [[8l-fT0l. 
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Nomenclature 

aij elements of covariance matrix of stochastic vectors a, and b. 

P^, Pk variance errors of the parameters before and after iterative update at 
iteration k 

Qk^ Q covariance error of the parameters at iteration k, constant covariance 

Rk, R^‘'' covariance error for measurements at iteration k, constant covariance at 

initial parameters i 

Xk vector of estimated state variables at iteration k 

H, Hi; matrix of derivatatives of /z, at iteration k 

Kyt nxN matrix represents Kalman gain 

Mix) mean square error between z and h 

V, Vk vectors of N errors of the measurements, at iteration k 

O sum of the squared residuals between z and h 

{ki elements of noise vectors w, and v 

^(t) conductivity of escaped ions at time t 

a, b two stochastic vectors 

A, B, R three constants 

D diffusion parameter 

h vector of N predicted measurements 

L side length of bone sample 

N, n, M number of measurements, state variables, iterations 
p a priori estimated relative variance 

T transpose 

t time 

w, Wk vectors of n perturbations of the state variables, at iteration k 
wd, wb noise of parameters D, and B 

X, X, xq vectors of n unknown, perturl^d, initial state variables 
z, Zk, Zav vectors with N measurements, at iteration k, avergare z 




Kalman filter can be used to determine material parameters from uncertain and 
inaccurate measurements. Aoki et al. [im used Kalman filter to identify Gurson’s 
model constant. They found that the accuracy of parameters prediction is affected 
by both specimen geometry and measurement type, and the shape of the tested 
specimen affects the convergence of the parameters. Also, they noticed that the 
rate of convergence can be improved by combining measurements of two differ¬ 
ent specimens in shape. The identification of Gurson-Tvergaard material model 
parameters via Kalman filtering technique is studied by Corigliano et al. ffT^ . 
They stated that the estimated values of the parameters are in well agreement with 
those obtained in previous work, but the initial suggested values for the seeking 
parameters affects the estimated parameters. 

Nakamura et al. [fT3ll implemented Kalman filter to determine elastic-plastic 
anisotropic parameters for thin materials using instrumented indentation. They 
observed that the initial chosen values for the parameters converged to a specific 
small area, but not to one point. Also, based on the convergence intensity, the 
parameters are determined. The same findings are obtained by using Kalman 
filter to determine the nonlinear properties of thermal sprayed ceramic coatings 
m. Bolzon et al. [fTSll used Kalman filter to identify parameters of a cohesive 
crack model. They reported that almost a linear correlation between convergent 
parameters is found, and the reason for the multiple local minimum might be 
related to using the linear Kalman filter for non-linear model. 

Vaddadi et al. [|T^ used Kalman filter to determine critical moisture diffusion 
parameters for a fiber reinforced composite. They estimated the parameters from 
the intensity of the convergence, which found to be in consistent with known 
values. Another study made by Vaddadi et al. IfTTlI to determine hygrothermal 
properties in fiber reinforced composite using Kalman filter. The parameters are 
extracted by reading the intensity of convergence plot. 

Kalman filter is an efficient way to filter noisy experimental data for determi¬ 
nation of material parameters. However, the initial suggested parameters required 
for Kalman filter should be chosen carefully, to avoid false local attractor. Also, 
the covariance error for the parameters noise almost assumed to be zero, which 
slow the rate of convergence and might lead to more than one intensity area for 
the predicted parameters. 

In this study, a methodology will be applied for using Kalman filter to deter¬ 
mine material parameters using uncertain measurements. The methodology starts 
by a way based on the mean square error to choose appropriate initial suggested 
parameters required for Kalman filter, and followed by a suggested way to choose 
the covariance errors for both state and measurements. The determination of dif- 
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fusion coefficients in bovine bone for generated data with different noises scatter 
from known parameters will be applied as a case study. A real measurements will 
be used also. 

2. Methods 

2.1. The Model 

Assume that an experiment resulted in N measurements obtained at differ¬ 
ent times, locations, temperatures etc. These are collected in a vector, z, with N 
measurements. The experimental data may be obtained at different known times, 
locations, temperatures etc. Measurements and all other data are available a pri¬ 
ori. 

In an attempt to predict the measurements a model, h = h{x) is used, with 
h being a vector of N predictions of observations. Further, v: is a vector of n 
unknown parameters defining the model based on variables such as position, tem¬ 
perature, time, etc. The unknown model parameters may describe the state of the 
system regarding, material, geometry or similar. In the present study, x is limited 
to parameters describing the material. 

Measurements always include systematic and non-systematic errors due to in¬ 
strumentation, indirect observations, gauges sensitive, irrelevant external influ¬ 
ence, and similar. Material parameter is sought but the experimental method may 
require a state parameters to be determined as well. Further, material parameters 
contain non-systematic errors due to thermal fluctuations, unstable structural con¬ 
figurations such as mobile dislocations, impurities, inclusions, unstable chemical 
composition, etc. Also inevitably, there is a difference between model and reality 
while a model never gives an exact description of the physical processes. Under 
ideal conditions the model would perfect in the sense that z = h{x). Here, only 
non-systematic errors or noise is considered. The model is defined for measure¬ 
ment i as 

z = h{x) + v, (1) 

where v is a vector with N errors due to inaccurate measurements z. The instant 
parameter x corresponding to the individual measurement i includes noise accord¬ 
ing to 

x=x -\-w, (2) 

where w is a vector with n errors caused by the parameter deviations. The elements 
of V and w are assumed to be uncorrelated. All elements of w and v are supposed 
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to be random, having the same respective stochastic distribution and for both a 
vanishing mean value is expected, cf. [|2l. 

Assuming that a set of parameters Xk is an estimate in the neighborhood of 
X, an improved estimate Xk+\ may be obtained by using linearized using a Taylor 
series which gives 

h{xk+\)-h{xk) + 'R{xk){Xo - Xk ), (3) 


when quadratic and higher order terms of x are neglected. On matrix form in¬ 
volved variables are 





r 

1 



h{x)= 


, H(x)= 

dxi 

dXn 

, X = 












dxi 

dXn 




Here, H is an A X n a Jacobian matrix. 


2.2. Least-Squares 

The system is supposed to be overdetermined, meaning that the number of 
measurements N exceeds the number of unknown parameters n. The best fit in the 
least-squares sense minimizes the sum of squared residuals with a residual being 
the difference between an observations z and the predictions h{x). The unknown 
X is obtained in a series of iterative improvements of the approximation x ~ Xk, 
where k is the iteration number. The initial parameters xq may be an educated 
guess based on previous measurements, data from resembling materials or other 
similar expectations. 

Solutions for non-linear systems (see Appendix A) may be obtained itera¬ 
tively. As an example, the Newton-Raphson method applied to these solutions 
gives the following recursive scheme, 

Xk+i =Xk + (H[Hi)“^Hj{z - h(xk)}, (5) 


where is an A x n matrix and a function of Xk- If convergence is reached, a 
local minimum of the sum of squared residuals has been found. To find global 
minimum, additional steps have to be taken. The drawback of the method is that 
convergence is not necessarily reached and is less likely if the a priori information 
of X is vague, unreliable or even misleading. This is especially accentuated when 
the measurements are noisy. Further, all data have to be present a priori. Modi¬ 
fications have been developed that allow an incremental treatment, which may be 
useful if data is continuously added (cf. IITSl f. 
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2.3. Kalman Filter 

The method of least-squares does not take the properties of the noises v and 
w as regards expeeted mean value and distribution into aeeount. Opposed to that, 
the Kalman filter is developed to use information about the noise a priori or as the 
measurements are assembled. The filter is an improvement of the least-squares 
method as it reeursively optimizes the unknown model parameters in seareh for 
least sum of squared errors. The method is developed with the partieular endeavor 
to effeetively handle noisy input data [l2]|. It is also operating inerementally so 
that new may be added as they appear in during ongoing measurements without 
loss of aeouraey. However, in the present study all measurements are supposed 
to be available when the optimization is initiated. The Kalman method is gener¬ 
ally assumed to be an effeetive method to filter noisy data eombined with a high 
eonvergenee rate. The derived algorithm is taking the eharaeter of the noise into 
eonsideration. The expeeted vanishing mean values for y and w are explieitly 
utilized. 

Assuming that the measured data is the predieted data based on the optimum 
material parameters with the addition of noise, as given by Eq. Q it is here 
assumed that the noise v a distribution with a zero mean value. 

Initially a set of parameters xq is seleeted based on a priori information from 
other measurements under same or similar eonditions or otherwise known data. 
From this, the parameters are iteratively updated using an algorithm on the same 
form as the least-square algorithm (of. Eq. as follows 

Xk+\ = Xk + - Kh)} for k=0,1,2, ...M , (6) 


where is an n x A matrix denoted the Kalman gain, and M is the number of 
iterations. In Eq. (|^. Normally the Kalman algorithm operates on single mea¬ 
surements one by one so that the Kalman gain is updated for every new measure¬ 
ments. This may be neoessary for interaetive proeesses where the aetion requires 
knowledge of the instantaneous state of the system. In the present study, iterations 
are simultaneously utilizing all measurements. The derivation of the Kalman gain 
(see Appendix B) is based on measurements added reeursively as in the original 
form of the filter. 

The optimal K^. that minimizes the mean square error is given by 
Ku = Pinl{nuPinl + Ru)-\ (7) 


where P^ is the varianee of the errors of the parameters before the iterative update, 
and is A X A matrix introduees the eovarianee error of the measurements that 
eomputed as 
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Ru = Var(v;t) ■ 


( 8 ) 


The is expressed as 


k+\ 


{I - ^k^k)Pk + Qk ’ 


(9) 


where Q^^ is an nxn matrix represents the eovarianee errors for the state parameters 
that is eomputed as 


0^ = 2Var(wi), 


( 10 ) 


The derivation of , Rk, and 0^ are given in details in Appendix B 


The algorithm involves reeursive use of the Eqs. Q, (|7]) and Q. The mea¬ 
surements eonsidered in eaeh reeursive eyele may be everything from a single 
measurement to all measurements. For non-linear problems, eaeh eyele may be 
repeated until eonvergent results are obtained. As an alternative, the entire reeur¬ 
sive seheme may be restarted and the resulting parameters P~, from previous 
applieation of the seheme are used as initial parameters. Also 0^ and Rk may be 
adjusted based on the improved information that is obtained. In the present study, 
all measurements are plaeed in a single set with N measurements, meaning that 
= z is a eonstant veetor with N elements. The number of reeursive eyeles is 
M and k = 1,2,..., M where M is preseribed or eonditional. The reeyeling is per¬ 
formed to aehieve a eonverged result for a non-linear problem. In eaeh eyele, Xk, 
H^:, Pk and therefore K^. is updated. 

By putting Rk = RI in Eq. Q and then taking the limiting result as 7? ^ 0 one 
obtains After inserting this into Eq. ([^ it is 

readily seen that the result is identieal to that of the non-linear least square method, 
of. Eq. ([^. The result is independent of and oonsequently also independent of 

Qk 


3. Methodology 

In [[II1[I1[I7]|, the Rk was ehosen to be a small peroentage of the measured 
data, and the 0^ value was ehosen to be zero. In a reeent study, the Rk value was 
ehosen as the differenee between the measured data and a predioted data a round 
the measured data, and the 0^ was ehosen to be unity llT9ll . 

To use the Kalman filter the parameters 0^ and Rk, and initial values for xq 
and Pq have to be defined. In the following, different strategies for ohoosing these 
values is desoribed. 
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For common usage of the Kalman filter, the ehoiee would he = Var(yjt) 
and Qi^ = 2Var(w^) aeeording to Eqs. ([^ and ( [T0| ). These variations are assumed 
to be known a priori, at least approximately. The information may be based on 
expeetation or derived from the present measurements, using a large variety of 
hypotheses. When the method is used reeursively, the indiees k allow for using 
Qj. or Rj, or both as funetions of time, position, ete. 

The Kalman gain, given by Eq. 0. with the selected Qj^ or R^, minimizes the 
squared error of the estimate Xk of x. A eondition for the derivation is that h(x) 
is a linear funetion of the parameters x. In the present study, h{x) is a non-linear 
function of x. The aim is to formulate a strategy for seleeting the free parameters 
Rk and not necessarily aeeording to Eqs. ([^ and ( [TO] ), so that the square of 
the error 

Mix) = ^YRr{z-h(x)), (11) 

is minimized. The study does not attempt to be exhaustive and the eonelusions 
are empirieal and based on a ease study. The seleetion of method is primarily 
based on eonvergenee rate. With a wide range of starting values, oeeasionally the 
estimate eonverges outside the range of interest and there is also the risk of failure 
in produeing eonverging results at all. These, disadvantages are also eonsidered 
in the seleetion of a suitable proeedure for seleeting R^ and 

In the first part of the study, the initial parameters xo are seleeted to eover a 
several orders of magnitude wide range of values. Under normal eireumstanees, 
this eannot be done for non-linear phenomena or realistie geometries or anything 
else for whieh an analytieal solution eannot be found, whieh may be the general 
ease. When the predictions are based on non-linear numerical calculations of field 
problems, e.g., using lengthy finite element analyses, usually only a few initial 
parameters xo ean be eonsidered. Here however, a wide range of initial parameters 
is examined as regards the mean square error At(vo). The aim is to obtain an 
overall pieture of the possibilities of fast eonvergenee or difficulties beeause of 
present loeal minima, sadle points ete. 

The seeond part is, the using of a suggested method and eomparing the re¬ 
sulting eonvergenee rate, the ability of produeing eonvergent results, and the per- 
eentage of eonvergent results from a range of the initial parameters with three 
additional methods. All four methods are deseribed in isubseetion 3.11 

3.J. Selecting Qi, and R 

In the present study, the results of constant 
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Qi^ = 2p^xx^ including Q^. = 0, 


( 12 ) 


are evaluated. The parameter p is an a priori estimated relative varianee. 

The relative varianee p of the state parameters is a quantity that possibly ean 
be guessed with more or less aeeuraey. Sinee the real state parameters jc are not 
known a priori eonstant may be either too big or too small depending on how 
aeeurate initial guess of jcq is. However, sinee the knowledge of the state variables 
inereases as the iterations proeeed, to stiek a eonstant may unneeessarily slow 
down the eonvergenee rate. As the iterations proeeed, the estimate Xk is improved 
and a better estimate for ean be used. Here, an updated Qk is also evaluated. 
With the variation of x being Wk = px, the unknown x is here assumed to best 
approximated with x^. The updated value given by Eq. ( [T0| ) then beeomes 

e, = 2/Var(x,). (13) 


The Rk taken as 

= V®/, i = 1,2, ....1681, (14) 


whieh is assumed to be a reasonable approximation of Eq. 0 . where measure¬ 
ment noise is given by 

= max|Var(zA, - %(4'^))l] (15) 


where max denoted to the maximum value, and i is the number of initial pa¬ 
rameters. 

In the suggested method, the Qk and are ehosen for large values based on 
the evaluation aeeording to Eqs. Q and Q, of. [seetion 5| for explanation. There 
are additional two Kalman filter methods, one method uses p = 0 in Eq. ( [T^ , and 
the another method uses p = 0.01 that gives the smallest standard deviation and 
largest peroentage of eonvergenee among different values for p, the mean values 
are elosest to eaeh other (see|Eig. 4). These two methods use Rk as 


R = [max|Var(zA? - hj^{x))\]I 


(16) 


so, the R is the largest squared element for the varianee of the differenee be¬ 
tween measured data, z, and predioted data with noise. Here, Q and R without k 
index denote eonstant eovarianees during iterations. 

The fourth method is the non-linear least square method, whieh is obtained by 
letting ^ 0 as it is deseribed in|subseetion 2.3 
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3.2. The Initial Predicted Parameters Error 


The initial predicted parameters error is an n x n matrix that contains the 
expected values for the errors between the unknown parameters xj, and the initial 


predicted parameters xq before the first iteration (see Appendix B). is chosen 
as 


Po = Var(x 




(17) 


where and are the two vectors that contain maximum and minimum 
of the selected initial values for the parameters. The choice for is the same in 
all cases in the present study. Only one large initial value is tested since the Pq 
introduces the variance between the seeking parameters and the initial parameters. 


3.3. Summary of Methods 


Table 1: Methods used 


Method 


Q 


Kalman filter 1 
(suggested method) 
Kalman filter 2 
Kalman filter 3 
Non-linear least square 


2 = Po as in Eq. (23 


Eq. (12 with p = 0 


Eq. (12 with p = 0.01 


R P 


Eq. ( 

26 

Eq. ( 

17 

Eq.i 


Eq.i 

17 

Eq.i 

16 

Eq.i 

17 


0 


D.N% 

92.09 

80.90 

80.19 

40.93 


4. Case Studies 

The determination of the diffusion coefficient of mammal bone using uncer¬ 
tain data is chosen to be a case study. Diffusion has recently been suggested to 
play an important role in transporting substances from the inner boundaries to the 
outer boundaries of bone. Therefore, knowing the diffusion coefficients in human 
bone are important to give required information for design of individual exercise 
programs that maximizes bone remodeling and bone healing. 

In the present study, the proposed four methods are applied to several simu¬ 
lated cases of generated data and the most effective method is applied to a case 
of real experimental data. The real experiment measures the amount of ions that 
leaves bovine bone samples that were put into a container with distilled water. 
During elapsing time the conductivity, ^(t), of the water increases in proportion 
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to the ionic concentration. The experiment is reported in Lindberg et al. [[T9l . 
Cubic bone samples with the side length L =10.1 mm from a bovine long bone 
were used. The concentration were measured using a SevenEasy S30 conductiv¬ 
ity meter from Mettler Toledo with an accuracy of +0.5% of the measured value. 
Further details regarding the experiment is found in [[T9l . 

The following model is suggested for the conductivity as a function of time, t. 


8 t 

m ^ , exp{-(2m - 1)VD-} , (18) 

m=\ '' 

where D is the diffusion constant, and A and B are unknown constants. The model 
is based on Pick’s law using relevant boundary conditions. 

The constants A and B provide the relation between the concentration in the 
bone sample and the conductivity in the distilled water. Putting t = 0 gives 

A = m + B, (19) 


where ^(0) taken to be the conductivity measured at t = 0. By using h(t) = 
^(t) - ^(0) as the measured quantity, the following model is obtained, 

M o ^ 

h(t) = B{l-J] ^2nm-iy ’ (20) 

m=l '' '' 


The two remaining constants D and B is determined using the Kalman filter. More 
details are found in Lindberg et al. IflPII . 

To the experiment is added a Monte Carlo set of 250 fictive measurements 
where generated data is used as measurements. The measured data is generated 
using the exact model Eqs. Q and (20) to compute the measurement vector z. To 
provide realistic conditions, a variation is added to the state parameters D and B 
and to the measurements. To this end, a random noise of 5%, i.e. \wd\ < 0.05D 
and \wb\ < 0.05B, is added. In the same way, noises v of 10%, 50% and 100% are 
included in the generated measurements z, i.e. z = h{t) + v, where v = Kh) 
and q equals 0.1, 0.5 and 1. Summation is performed for the 24 different times 
of measurement. Further, the wd, wb and v are uncorrelated and the probability 
density is constant within the limits of the respective noise. Also the noises for 
each measurements are uncorrelated. Thus, 


Co\{§i, '&j) 



( 21 ) 
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where and represents the elements of eaeh noise veetor w^, Wb and v at 
individual measurements i and j. 


As observed from Eq. (20), the only available time unit is provided by L}/D. 
A typieal experiment lasts for around 11.7 time units, i.e., 0 < t < n.lL^ID with 
generated measurements taken in time intervals of around 0.0032L^/D. 

In an attempt to obtain less non-linear formulation of Eq. ([20]), a semi-linear 


model that is obtained after simplifieation and taking log of both sides of Eq. (20) 
one obtains 


h*{t) = B* + f{D,t), 
where 

h*(t) = log m)}, B*=log(B), 
and 

M g 


( 22 ) 


f{D,t) = log 




exp{-(2m - ifn'^D—] 


^ 7T^{2m - 1 )' 


whieh makes h*(t), a linear funetion of B* while f{D, t) is known to be rather 
small. Note that log(l - e~^) e~^ as x ^ oo and therefore f{D,t) deeays 

exponentially for large with inereasing time whieh makes the linearization with 
respeet to D undoable. 


5. Results and Discussions 


The aeeuraey of the series in Eq. (20) was studied in this work to reduee the 


eomputational time, it showed that the number of terms m may be ehosen to be 
around 200 terms to obtain an aeeuraey of 99.9% in the middle of the sample at 
t = 0. Already at the seeond measurement at t = 0.003L^/D only four terms are 
required to obtain the same aeeuraey. Eor systems with large amounts of data, a 
strategy for the seleetion of the number of terms in the series eould save eonsid- 
erable eomputation time. However, here this is not the primary foeus and the the 
ealeulations where not very mueh time eonsuming, whieh is why all ealeulations 
where made using M = 200 terms. 

Measurements from known D and B using wd = O.ID, wb = 0.15, and 
V = O.lzov where Zav is the average of all generated measurements i.e Zav = Zi- 
Then the Kalman filer was used to obtain approximations of D and B. This was 
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done first for the non-linear model Eq. ( [20| ), and then for the partly linearized 
model Eq. (22). To eompare the aeeuraey in finding the least square error of the 
predieted measurements based on dilferent state parameters found using model 
Eq. ( [20| ) and model Eq. ( [22] ) respeetively, the least square error M{x) was ealeu- 
lated using Eq. The non-linear model gave A\{x) = 15% while the partly 
linearized model gave 30%. Eor the reason of this, the partly linearized model is 
given up and the model Eq. (20) is used in the eontinued analysis. 

Assuming that the information of the parameters D and B is uneertain, the 
Kalman filtering has to eonverge from initial values that are several orders of 
magnitudes dilferent from the true values. To explore what this means, the mean 
square error is ealeulated for a very wide range of initial values. Mean square 
error here refers to the error of ealeulated estimates of the initial guesses Dq and 
Bq direetly eompared with the true values of D and B. No iterations are made. 
To eover a large variety of initial values a mesh 101 values for D and for eaeh 
101 values for B are used. The mesh eovers values of Dq ranging from 0.00ID to 
lOOD and values of Bq from 0.015 to 1005. [Eig. Ij shows the mean square error 
for the initial values of Dq and Bq for the eonstrueted mesh. 



Fig. 1. The mean square error of the eonduetivity predietions for dilferent initial 
values of the relative initial parameter values Dq/D and 5o/5. The noises are in 
the range +5% for parameters and +10% for measurements. 
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Regions with small gradients are found along a line for whieh DqID approxi¬ 
mately equals The several loeal minima along this trajeetory eorrespond 

to the resolution of the grid and are merely graphieal anomalities. Also in the re¬ 
gion of small Dq and small Bq the gradient is very small. This makes the eonver- 
genee rate of any gradient driven algorithm small. However, no loeal minimum 
seems to be present and eonvergenee should be possible in entire range of initial 
values even if the eonvergenee may be very slow in the above deseribe regions. 
A elear minimum mean square error is found around the elose to D and B. The 
result is strongly influeneed by the rapid ehanges due to the exponential behavior 
ofEq. ( [20] ). 

For the study of the Kalman filter measurements z are generated for 250 mea¬ 
surements in the time interval 0 < t < ll.lL^/D. In this study, D and B are known 
parameters, but in a real ease they are not known but believed to be in the neigh¬ 
borhood of the a priori guess. As it is suggested by |Fig. 1[ a very large variation 
of eonvergenee rates are antieipated. 

The Kalman filter is studied by gathering 10 generated experiments together, 
eaeh experiment has 25 measurements, and eaeh measurement has random noises 
of +5% for parameters and +10%, ±50%, and +100% for measurements, whieh 
eonstruets a measurement veetor Zi, i=l,2.3.,250. The initial seleeted state vari¬ 

ables are ehosen from O.IZ) to lOD for Dq and from 0.15 to 105 for Bq. 

The initial parameters error Pq is seleeted aeeording to Fq. ([T7|) as follows 


9.9^Dl 0 

0 9.9252 _ ’ 


The is updated using Fq. (|^. 


(23) 


5.1. Different R and Qt 

First values are ealeulated as in Eq. ( [T^ . So, is the squared largest 
element of the varianee eneountered so far, i.e. so far means before iteration k. 
The effeet of Q^. is studied by using. 


(Dk - Dof 0 
0 (5, - 5o)2 


(24) 


This means that inereases as the iterations proeeed. The expeeted effeet is that 
the eonvergenee rate inereases with inereasing Whieh is seleeted to speed up 
the eonvergenee rate. 
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Fig. 2 a) shows a color plot of the obtained values for Di/D, here 1 is num¬ 


ber of iterations. The markers (x) that are included show the obtained values of 
Di/D and Bi/B using Kalman filter, for 41x41 initial values Dq/D and Bq/B. The 
obtained {DilD,BilB) outside the seleeted range are exeluded from eonvergenee 
plot to make it in the same range as the eolor plot. The white areas in the fig¬ 
ure give 21.5% of the obtained Di/D that found to have negative values with no 
physieal meaning. Consequently, this leads to divergent result in following itera¬ 
tion since the exponential term in the diffusion model would have large positive 
values. A large step for Kalman gain is seems to be the reason, whieh ean 
be forced to be small by assuming R large, but the rate of eonvergenee would be 
slow. 



d„;d d„/d 

Fig. 2. Color plot of Dk as a funetion of Dq and Bq. The markers (x) show the 
position of the resulting Dk and Bk for the 1681 initial starting points in the 
range, (a) k = 1 iteration, and (b) k = 20 iterations. Noises are +5% for the 
parameters and +10% for the measurements. 


Fig. 2 b) shows the eolor and eonvergenee plot after 20 iterations. It ean be 


seen that using larger values for Qj. inereases the rate of eonvergenee but the di¬ 
vergent eases inereased from 21.5% to 22.4%. 

The effeet of R values on the eonvergent T>i/T> is shown in |Fig. 3[ a). The 
figure shows that the percentage of the number of eonvergent D\ID over the total 
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number of D\ID almost around 80% as the R value increases from 0.0000395^ 
to 0.6B^, after that it begins to increase to 95% as the R increases to 3B^, then it 
almost keep around this percentage with larger values for R. A possible reason is 
that the Kalman gain step becomes small as the R value increases. 
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Fig. 3. Effect of R on Di/D. (a) convergent of DijD over total number of DijD. 
(b) standard deviation and mean value, the standard deviation values over 5 are 
excluded from the graph. Q = 0. 


The effect of R on standard deviations and mean values are shown in Fig. 3 'b). 
The figure shows that the standard deviation decreases as the R increases untill 
R = 1.55^, then it increases slowly, also the mean value decreases as the R in¬ 
creases to 0.025^, and then increases again. The R/B^ (Kalman) in the graph 
denoted suggested R for ordinary Kalman filter based on the average value of z 
and the 10% measurements noise, i.e R = {Q.lzavf, which is not the optimal case. 
One reason is that standard deviations and mean values are computed after one 
iteration, also the negative values of D are included. Slightly different results is 
obtained after more iterations, reason should be the non-linear model. 

Fig. 4(a) shows the effect of different values of Q using 0.01 < p < 1 and 


R as given in Eq. ( [T^ on the convergent numbers of D^jD over total number of 
D^/D. It is clearly seen that the number of convergent D^/D decreases as the p 
increases. Using large Q with small R affects the Kalman gain stability, this could 
be a possible reason. 
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Standard deviation of D^D 
Mean value of D /D 


p=1, Q=2var(w) 
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Fig. 4. The effect of p on D^jD. (a) on convergent of D^ID over total number of 


Dj/D (b) standard deviation and mean value of D^/D. R is chosen as in Eq. (161. 


Fig. 4 b) shows the effect of p on the standard deviations and mean values 


for the obtained D^jD. The figure shows that standard deviation increases as p 
increases, and mean values play around 1.7 as the p increases. One reason could 
be that using large values of Q with small values for R increases the Kalman gain 
step. The Q value at p = 0.71, random walk (cf. [|2l), and at p = 1 , the suggested 
variance as in Eq. ( [T3| ), does not give significant effect on both standard deviation 
and mean value. A non-linear model might be a reason, so it is expected to have 
effect with more iterations. The p value is chosen for 0.01 that give smallest 
standard deviation and largest percentage of number of convergent D over total 


number of D to be one method (see Table 1). 


5.2. Large R and Q 

By studying Eqs Q and Q, and |Fig. 3| it is believed that choosing R large, 
decreases the risk of ending up with divergent results for Dj,. |Fig. 5| ^a) shows that 
increasing p for large values using R large increases the percentage of the number 
of convergent Dj by 3%, this is an indication that using p > 0 with R large covers 
a wide range of initial parameters that converges. 
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Fig. 5. The effect of p on D^jD using R large, (a) convergent of D^/D over total 
number of D^/D (b) standard deviation and mean value of D 5 /D. The large R is 
chosen as in Eq. (|2^. 


The effect of Q on the standard deviations and mean values of the obtained 


D 5 /D are shown in Fig. 5 'b). ft can be seen that the standard deviation and mean 


value decrease as the p increases, untill optimal is obtained with p > 7.8. Also, 
increasing p speeds the rate of convergence. Again, this could be the effect of 
non-linearity of the model. 

ft can be concluded that using R and Q large increase the possibility of conver¬ 
gence for a wide range of Dq and Bq and speeds the convergence of parameters. To 
choose R for large value, a vector contains the largest variance value between 
measurement, z, and predicted , h{xo), for each initial parameters is chosen as 




(25) 


where s is the number of initial parameters, i.e s = 1,2,...., 1681, and values 


are given by Eq. (15), then the R value is chosen as the maximum value in the 
whole initial parameters combinations as 


R = max(Vi)/, 


(26) 


To choose Q large, the Q is chosen to be equal to Pq as in Eq. (231, since the 
Q represents the variance between the seeking parameters and instant parameters. 
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Fig. 6. Color plot for Dk merged with eonvergenee plot between {Dk,Bk) for 
initial starting points Doand 5o after (a) 1 iteration, for large R and Q = 0, and 
(b) 20 iterations for large R and large Q. Q is seleeted to be large and equal to 
as in Eq. (231, while R is ehosen as in Eq. (261. Noises are +5% for the 
parameters and +10% for the measurements. 


5.3. Different methods 

The method with the seleetions of Q and R, is eompared with with the seeond, 


the third, and the fourth method that are explained in subseetion subseetion 3.1 


and summarized in[TJ The R in the seeond and third methods is ehosen as the max¬ 
imum differenee between generated data, z, and generated data with the noises, 
h(x). The Pq for the seeond and third method is ehosen as in Eq. (231, the same as 
in the suggested method. Eig. 7 'a) shows that 92.1% of the initial seleeted Dq/D 
and Bq/B give eonvergent results for diffusion eonstants. Around 98.3% of those 
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values are eonverged to 1.042 after 50 iterations, while the 1.7% left (blue area 
in the eolor plot) almost eonverged to the same point after 180 iterations. The 
Kalman filter with p = 0.01 gives almost the same eonvergent area of D^q/D as 
the Kalman filter with p = 0 around 80%, but it speeds the rate of eonvergenee as 
shown in |Fig. 7( e), while Kalman filter with p = 0 needs a large number of iter¬ 
ation for the parameters to eonverge as shown in |Fig. 7[ b). The non-linear least 
square speeds the rate of eonvergenee but with 40.93% of the initial seleeted pa¬ 
rameters that give eonvergent results for diffusion eonstants as shown in|Fig. 7|^d). 
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D„/D d„;d 

Fig. 7. Color plot for D^^jD merged with eonvergenee plot between 
{D 5 q/D,B 5 qIB) for initial starting points Dq/D and Bq/B after 50 iterations for: 

(a) suggested Kalman filter with large R and large Q, (b) Kalman filter with 
p = 0, (c) Kalman filter with p = 0.01 , and (d) non-linear least square method. 
For generated data with noises of +5% for parameters and +10% for 
measurements. 

The suggested Kalman filter eovers a wide range of the initial seleeted values 
for Do and Bq that give eonvergent results for diffusion D 50 eompared with Kalman 
filter with p = 0, Kalman filter with p = 0.01, and non-linear least square. This 
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makes it an appropriate method to determine diffusion coeffieient if the a priori 
information is rare or there is a large variation in the parameters, such as bone as 
an inhomogeneous material. Also, it speeds the rate of convergence. The Kalman 
filter with p = 0.01 might be a good choice if the a priori information in the 
adjacent area of the seeking parameters are enough. Also, the non-linear least 
square might be one choice but for small range of initial selected values for Dq 
and Bq. 

The effect of different ranges for Dq and Bq against the standard deviations 


and mean values for the four methods after 5 iterations are shown in Fig. 8 The 
ranges are selected to start from a small range and end with a large range as 


1 

Tb 


< selected initial ranges < iD, 


iB 


< selected initial ranges < iB, 


where i = 1.2,1.6,2,., 10. The figure shows that the standard deviations and 

mean values increase as i increases. A reason for that could be that the number of 
initial parameters increases as the ranges increase. It can be noticed that Kalman 
filter with p = 0.01 gives the smallest standard deviations while the suggested 
method gives standard deviations closest to Kalman filter with p = 0. On the 
other hand, the mean values obtained by suggested method is found to be the 
smallest among the other three methods. This is an indication that the suggested 
method and Kalman filter with p = 0.01 might converge around the same speed. 
On the contrary, the non-linear least square gives very large values for standard 
deviations for i > 1.6, these values are above 2 and excluded from the graph as 
shown in Fig. 8 'a), and gives large mean values for D^/D as shown in Fig. 8[ b). 
A possible reason for that is the divergent result that obtained after 5 iterations, it 
has some large negative values for D^jD. 
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Fig. 8. Standard deviations and mean values for obtained D^/D against ranges of 
seleeted initial parameters for the four methods after 5 iterations, (a) standard 
deviation and (b) mean value. For generated data with noises of +5% for 
parameters and +10% for measurements. 


Fig. 9 


shows the standard deviations and mean values of the obtained Dk, 
k = 1,2,...., 100 against number of iterations for large initial range, i = 10, for 
the four methods. The standard deviation values for for the suggested Kalman 
filter deerease slowly up to 70 iterations, and the mean value inereases as number 
of iterations inereases to 70 iterations, and almost stable after that. The reason 
for that is some of the obtained values of D/t were stuek to small values up to the 
70 iterations, moving with small steps of Kk until they reaehed to good predieted 
values to move with larger steps. 
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Fig. 9. Standard deviations and mean values for obtained against number of 
iterations for the four methods after 5 iterations, (a) standard deviation and (b) 
mean value. For generated data with noises of +5% for parameters and +10% for 
measurements. ( points over 2 and less than 0 are exeluded from the graph). 


The standard deviation and mean value using the Kalman filter with p = 0 
deerease slowly as the number of iterations inereases. This is expeeted sinee R in 
the seeond method was ehosen to a small value and Q vanishes, whieh makes the 
Kalman gain is small. On the other hand, the standard deviation and mean values 
of Dt for Kalman filter with p = 0.01 and for non-linear least square eonverge 
quieker than the others. This is an indieation that methods 3 and 4 ean be used 
effeetivly if the apriori information fortunately adjaeent to the seeking parameters. 

The determination of diffusion eonstants for generated data with ±5% param¬ 
eters noise and ±50%, and +100% measurements noise almost follow the same 
trend as in +5% parameters noise and +10% measurements noise. [Table 2| shows 
the standard deviation, mean, and the pereentage of numbers of eonvergent D^olD 
over total number of D^o/D (D.N%) for the seleeted range with i = 10, for the four 
methods. 
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Table 2: Standard deviations, mean values and D.N% for the obtained D 50 /D for 
the four methods with noises of +5% for parameters and +10%, ±50%, and 
+ 100% for measurements. 


Random noises 

Method 

Standard deviation 

Mean 

D.N% 

±10% 

Suggested Kalman filter 

0.09 

1.03 

92.09 


Kalman filter, p = 0 

0.49 

1.22 

80.90 


Kalman filter, p = 0.01 

0.03 

1.04 

80.19 


Non-linear least square 

8.71x10-'” 

1.04 

40.93 

±50% 

Suggested Kalman filter 

0.06 

0.93 

91.96 


Kalman filter, p = 0 

0.54 

1.06 

80.37 


Kalman filter, p = 0.01 

0.003 

0.94 

80.67 


Non-linear least square 

6.19x10'*® 

0.94 

45.63 

±100% 

Suggested Kalman filter 

0.03 

1.69 

92.03 


Kalman filter, p = 0 

0.66 

1.44 

85.24 


Kalman filter, p = 0.01 

0.22 

1.66 

85.66 


Non-linear least square 

0.064 

1.69 

42.47 


The table shows that the suggested Kalman filter gives a compromise results 
compared with the other methods. The most important difference is the percent¬ 
age of obtained D^q/D that give convergent results (around 92%), that is found 
to be large compared with the others, which means high possibilities to deter¬ 
mine the diffusion coefficients from the selected initial range even it was a large 
range, and for measurements noise up to +100%. The table also shows that the 
suggested Kalman filter speeds the rate of convergence as well as Kalman filter 
with p = 0.01 and non-linear least square but with higher D.N%. By comparing 
standard deviations and mean values as in the table, the standard deviations for 
Kalman filter with p = 0 give the largest values among the other methods, and 
the mean values are pretty a way from the convergent parameters. A reason for 
that is Kalman filter with p = 0 converges to many different points with small 
Kalamn gain step, which means more iterations are needed for convergence. In 
sum, the suggested method can be applied effectively for both rare and sufficient 
information about the seeking parameters, and can be applied for a wide range of 
initial parameters. 

The behavior of mean square error for a selected real bovine bone sample is 
found almost the same as the generated one in [Fig. 1[ The initial predicted param¬ 
eters Dq and Bq are chosen to be 41 x 41 combinations, with 0.001 < Dq < 0.1 
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mm^/min and 6 < Bq < 600 [iS/mm, which constructs a large combinations be¬ 
tween the initial parameters. The initial predieted parameters error Pq, eovarianee 
error for the parameters Q, and eovarianee error for the measurements R are eho- 
sen based on both suggested method and these eombinations. |Fig. 10| shows the 
eonvergenee plot for the two parameters D and B that obtained using the suggested 
Kalman filter for the real sample after 1, 10, 50, and 250 iterations. 



Fig. 10. Convergenee between and Bj^ obtained using suggested method for a 
real sample for 1, 10, 50, and 250 iterations 


The initial predieted parameters that give eonvergent results for Dj, are found 
to be 92.2 % from the eombination, that eonverged around a speeifie line after 
some iterations as in 10 and 50 iterations (see Fig. 10[ ). Around 92.5% of them 
eonverged to the seeking unknown parameters D = 0.0144 mm^/min and B = 
73.602 p,S/mm after 50 iterations, while the 7.5% left eonverged to the same plaee 
almost after 250 iterations. Small Kalman gain steps at the 7.5% might be a 
possible reason for the slowly eonvergenee. The eonvergenee line is expeeted 
sinee the mean square error has smaller values along this line as shown in [Fig. 1 


The eonduetivity versus time for analytieal model Eq. ( [20] ) using D and B ob¬ 
tained by the suggested Kalman filter and experimental data for the real sample is 
shown in |Fig. 11] The figure shows that the analytieal funetion fits very well with 
the experimental data aeeompanied with mean square error of 0.85 ([iS/mm)^. 
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Fig. 11. Conductivity versus time for analytieal and experimental for the real 
sample 


6. Conclusions 

Four methods to ehoose appropriate initial parameters, eovarianee errors for 
parameters, and eovarianee errors for measurements, required for Kalman filter 
for determination of material parameters, is investigated in this work. The meth¬ 
ods are applied to generated data with +5% parameters noise and +10%, ±50%, 
and +100% measurements noise for known parameters. The suggested method 
eovers a wider range of the initial suggested values for the parameters than the 
standard Kalman filters, and non-linear least square, whieh enhanees the possibil¬ 
ities of eonvergenee around the seeking parameters. The suggested method speeds 
the rate of eonvergenee eompared with the other methods. Very good results are 
obtained for diffusion eoeffieient in bovine bone as a ease study. 

Appendix A. Non-linear least squares 

The sum of the squared residuals is written 
O = (z --/?) (A-1) 


28 





By putting the derivative of O with respeet to x to zero, a set of n equations is 
obtained, 


50 

dx 


-=0 


H(A:r {z - h(x)} = 0. 


Insertion of Eq. Q into Eq. (|A-2|) gives the following. 


H^lz - h(M * rfnix - Xk). 


(A-2) 


(A-3) 


Provided that the n x n matrix H^H is non-singular, one obtains the following 
approximation, 

x~Xk + (HTH)-'H^{z - h(xk)}. (A-4) 


Eor linear systems the Jaeobian H is independent of x, whieh makes the solution 
in Eq. (|A-4|) exaet. 


Appendix B. Kalman filter derivation 

The varianee of the errors of the parameters before, P^, and after, Pj„ the 
iterative update are 

P^ =Yar{xk-Xk-i), and PA:=Var(x^ - . (B-1) 


The funetion Var(a) is an Z x / matrix with the elements aij given by the veetor 
a = (ai, as follows (of. [I20ll ) 

aij = E(XiXj) - E(Xi)E(xj ), (B-2) 


where is the statistioal mean value of the stoehastio variable &. Henoe, the 
matrioes P^ and are symmetrio with the dimension nx n. 

Using Eqs. (B-1) and ([^, a relation between P^ and Pt is obtained as. 


= Yar(xk+i - Xk) 

= Var(ji; + Wk+i- h) 

= VarCxjt -Wk + Wk+i- h) 

= VarCxjt - Xk) + Var(-w,t -h wu+i) 
= VarCxjt - Xk) + 2Var(wi) 

= Pk + Qk 
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where = 2Var(wi;) is an n x n matrix. Note that the eovarianee of the 
supposedly uneorrelated quantities {x - Xk) and Wk respeetively vanishes. 

On the other hand, substituting Xk in Eq. ( B-l[ ) by using Eq. (|^, Pk ean be 
expressed as 


Pk = Var[xi - {xk-i + ^k{zk - h{xk-i)})] . 

Replaeing Zk aeeording to Eq. Q leads to 

Pk = Var(X;t - 4-1 - ^k{h{Xk) + - /z(4-i)}) 

= Var(X;t - Xk-i) - KkCov(h(xk) - /z(4-i), Xk - 4-i) 

- CoY(xk - Xk-\,h{xk) - /i(4-i))K^ 

+ KkWds{h{xk) - h{xk-i))Kl 

+ K^Var(vOK[, 

(B-4) 

where the funetion Cov(a, b) gives the eovarianee of the stoehastie veetors a and 
b. The funetion Cov(a, Z?) is an m x m matrix with the elements aij given by the 
stoehastie variables a = {a^, ...a„,)^ and b = (bi, ...b,„y as follows (of. [l20ll ). 

aij = E{aibj) - E(a,•)£((?;•). (B-5) 

It is used that the eovarianee between the elements of h(xi) and Vk for any i 
and k vanishes. The Taylor series in Eq. Q giving h(xk) - h(xk-i) = Hkixk - 4-i) 
results in 


Pk = Pk - KkHkPk - P-kHjKl + KkHkPktilKl + KkRkKj 
= (/ - KkHk)Pk(^ ~ + KkRkKj , 

where Rk is the eovarianee error of the measurements and an A^x matrix defined 
as follows 

Rk = Yax(vk). (B-7) 


Sinoe the Xk and Vk are mutually uneorrelated, Pk beeomes a diagonal ma¬ 
trix that oontains errors between the parameters before and after an iteration, the 
Kk that minimizes the error ean be obtained by taking the derivative of the traee 
Tr(i*<;) with respeet to K^; and putting it equal to zero. Taking the traee of the first 
equality in Eq. (|B-6[) provides 
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Tr(n) = Tr(Pp - Tr(K,H,P,-) - Tr(P,-H^K^) + Tr(K,H,P,-H^K^) + Tv{K,R,Kj) 

(B-8) 


Using following identities for the matrix trace (cf. 11211 ) 


Tr(PrHX) = Tr(K,H,P:), 


ffTr{AC) 


=C^ and 


dTv(AFA^) 


=2AF. 


(B-9) 


dA ■ “ dA 

The last equality requires that the matrix F is symmetric. The derivative of Eq. 
( |B-8[ ) with respect to can be written 

dTv{P,) 


dKt 


= -2P;nJ + Hi + 2K,R 




k^k 


(B-10) 


The optimal Ki, can be obtained by putting the right term in Eq. [B-101 to zero 
as next 


K, = PrHl(H,PrHl+ /?,)-'. 


(B-ll) 


The Pk associated to the optimal K^. can be written as in Eq. [ B-12[ by substi¬ 
tuting Eq. [B-111 into Eq. ( B-6| ) as next 

Pk = {I- K,H,)F,- + {-P-,nl + K,(H,P,-Hl + Rk)]Kl 


Using Eq. ( |B-11| ) readily gives 
Pk = (i-K,n,)p-,, 
and consequently 

p;,,=(i-K,llk)P-k+Qk^ 


(B-12) 


(B-13) 
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